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ABSTRACT 

While the high-redshift quasar luminosity function closely parallels the hierarchical growth of dark 
matter halos, at lower redshifts quasars exhibit an anti- hierarchical turn-off, which moves from the 
most luminous objects to the faintest. We explore the idea that this may arise from self-regulating 
feedback, caused by quasar outflows. Using a detailed hydrodynamic simulation we calculate the lumi- 
nosity function of quasars down to a redshift of z = 1 in a large, cosmologically representative volume. 
Outflows are included explicitly by tracking halo mergers and driving shocks into the surrounding in- 
tergalactic medium, with an energy output equal to a fixed 5% fraction of the bolometric luminosity. 
Our results are in excellent agreement with measurements of the spatial distribution of quasars on 
both small and large scales, and we detect an intriguing excess of galaxy-quasar pairs at very short 
separations. Our results also reproduce an anti- hierarchical turnoff in the quasar luminosity function, 
however, this falls short of that observed as well as that predicted by analogous semi-analytic models. 
The difference can be traced to the treatment of heating of gas within galaxies and the presence of 
in-shock cooling. Calculations of the mass fraction of gas above the critical entropy show a strong 
redshift dependence with close to 20% of the baryons being above this limit at z = 1. Volume fractions 
show an even stronger trend with redshift and some sensitivity to resolution due to the tendency of 
high entropy gas to occupy low density regions. The simulated galaxy cluster Lx — T relationship is 
close to that observed for z m 1 clusters, but the simulated galaxy groups at z = 1 are significantly 
perturbed by quasar outfiows, suggesting that measurements of X-ray emission in high-redshift groups 
could well be a "smoking gun" for the AGN heating hypothesis. 
Subject headings: 



1. INTRODUCTION 

In the low-redshift universe, active galactic nuclei 
(AGN) are not very active. While at high redshifts, 
the quasar luminosity function increases with time, since 
z K, 2 the number density of optically-selected AGN 
has been dropping dramatically (Schmidt & Green 1983; 
Boyle et al. 1988 Koo & Kron 1988; Pei 1995; Boyle et 
al. 2000; Fan et al. 2001). Deep X-ray surveys have 
shown that this downturn occurs anti-hierarchically, such 
that the spatial density of AGN with higher X-ray lumi- 
nosities peaks earlier than that of lower-luminosity AGN 
(Stcffen et aL 2003; Ueda et aZ. 2003). Complementary 
emission line studies suggest that this trend is driven by 
a decrease in the characteristic mass of actively grow- 
ing black holes (Heckman et aZ. 2004), and is likely to 
closely parallel the formation history of early-type galax- 
ies {e.g. Granato et aZ. 2001, 2004). Furthermore, optical 
and near-infrared observations indicate that the largest 
galaxies were already in place by z sa 2, while smaller ones 
continued to form stars at much lower redshifts (Pozzetti 
et al. 2003; Fontana et al. 2004; Glazebrook et al. 2004; 
van Dokkum et al. 2004; Treu et al. 2005), and a sim- 
ilar trend is observed in the morphological evolution of 
galaxies (Bundy, Ellis, & Conselice 2005). 

Yet, despite these observations, such widespread 
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galaxy "downsizing" (Cowie et al. 1996) was unex- 
pected. The A Cold Dark Matter (ACDM) model, 
while in spectacular agreement with observations {e.g. 
Spergel et a/. 2003), is a hierarchical theory, in which 
gravitationally-bound structures grow by accretion and 
merging. Superimposed on this distribution is the bary- 
onic component, which falls into the dark-matter poten- 
tial wells, shock heats, and must radiate this energy away 
before forming stars (Rees & Ostriker 1977; Silk 1977). 
The larger the structure, the longer it takes to cool, and 
thus galaxy evolution should be even more hierarchical 
than structure formation. 

Recently, several theoretical studies have shown that 
the missing element in this picture could be kinetic feed- 
back from AGN. As bulge and black hole masses are 
closely related (Gebhardt et al. 2000; Merritt & Ferrarese 
2001; see also King 2003, 2005), such outflows would 
have the largest impact on the largest forming elliptical 
galaxies, suppressing their formation first (Scannapieco 
& Oh 2004, hereafter SO04; Binney 2004; Di Matteo et 
al. 2005; Croton et al. 2006). This would also help to ex- 
plain the X-ray luminosity-temperature relationship ob- 
served in the intracluster medium (ICM) in galaxy clus- 
ters. If nongravitational heating were unimportant, the 
gas density distribution would be self-similar, resulting in 
Lx oc (Kaiser 1986), but instead the observed slope 
steepens considerably for low-temperature clusters {e.g. 
David et al. 1993; Arnaud & Evrard 1999; Helsdon & 
Ponman 2000). Furthermore the 100 keV cm^ level of 
preheating necessary to explain this discrepancy (Cava- 
here et al. 1998; Kravtsov & Yepes 2000; Wu et al. 2000; 
Babul et al. 2002) is not arbitrary. Rather it corresponds 
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to the threshold value for cooling within the age of the 
universe (Voit & Bryan 2001; Oh & Benson 2003). 

Taken together, these observations are strongly sug- 
gestive of a model in which the properties of the ICM, 
the formation history of elliptical galaxies, and the evo- 
lution of the quasar luminosity function are all set by 
self-regulating AGN feedback. In fact, S004 have shown 
that the addition of AGN outflows into the semi-analytic 
model developed by Wyithe & Loeb (2002; 2003) can re- 
produce the drop in the AGN luminosity at low redshifts, 
as heating this gas slows the accretion of matter onto 
the supermassivc black hole. The central AGN engine 
is then starved of fuel and a strong suppression at the 
bright end of the luminosity function occurs. The main 
drawback of the model, however, is that precise match- 
ing of observational results requires fine tuning of the 
model parameters, such as black hole mass and outflow 
efficiency. 

Several other recent numerical and analytic investiga- 
tions have sharpened our understanding of various as- 
pects of this process. Binney (2004) conducted an ana- 
lytic study of the impact of AGN on inhibiting gas cool- 
ing in large galaxies. Di Matteo et al. (2005) carried out 
Smoothed Particle Hydrodynamic (SPH) simulations of 
AGN outflows in individual galaxy mergers and studied 
the role of feedback in determining the colors of elliptical 
galaxies and establishing the relationship between stel- 
lar velocity dispersion and black hole mass. This suite 
of simulations was later related analytically to the global 
evolution of quasars and elliptical galaxies in series of 
papers by Hopkins et al. (2005a, 2005b, 2005c, 2006). 
Levine & Gnedin (2005) combined cosmological simula- 
tions with an analytic model to constrain the filling factor 
of AGN outflows as a function of redshift. Scannapieco, 
Silk, & Bouwens (2005) emphasized the role that quasars 
may play in the downsizing of the star-forming galaxy 
population. Levine & Gnedin (2006) studied the impact 
of AGN outflows on the matter power spectrum. Menci 
et al. (2006) used a semianalytical model to study the 
role of AGN feedback on the color distribution of galax- 
ies from z = to z = 4. The importance of the nature of 
gas accretion in determining the effectiveness of feedback 
processes has recently been discussed in Dekel & Birn- 
boim (2006) and Cattaneo et al. (2006). Finally, Cro- 
ton et al. (2006) combined semi-analytic models with the 
dark-matter evolution taken from the Millennium Sim- 
ulation (Springel et al. 2005) to study the impact of a 
more temporally-extended model of AGN feedback on 
the bright end of the galaxy huiiinosity function. 

In this paper we undertake the first detailed hydrody- 
namic simulations of quasar outflows in a general cosmo- 
logical context. Adopting a burst model that associates 
AGN with merger events and a global outfiow model that 
is based on our simulations of high-redshift starbursts 
(Scannapieco, Thacker, & Davis 2001, hereafter STDOl), 
we are able track in detail the impact that quasars have 
both on their own formation, and on the properties of 
galaxy clusters and the intergalactic medium (IGM). 

The structure of this work is as follows. In §2 we de- 
scribe our overall numerical approach, method for quasar 
identification, and implementation of outflows. In §3 we 
compare our simulation results with measures of quasar 
clustering and the observed luminosity function in the 
optical and the X-ray bands. In §4 we study AGN feed- 



back in a more global context, examining its impact on 
the properties of the intergalactic and intracluster media. 
In §5 we present a discussion and conclusion. 

2. SIMULATIONS OF QUASAR FORMATION 

2.1. Overall Numerical Model 

Motivated by measurements of the cosmic microwave 
background, the number abundance of galaxy clus- 
ters, and high-redshift supernova distance estimates (e.g. 
Spergel et al. 2003; Vianna & Liddle 1996: Riess et 
al. 1998; Perlmutter et al. 1999), we focus our attention 
on a Cold Dark Matter cosmological model with param- 
eters h = 0.7, flo = 0.3, riA = 0.7, rib = 0.046, erg = 0.9, 
and n = 1, where h is the Hubble constant in units of 
100 km s~^ Mpc~^, ^^a, and fib are the total mat- 
ter, vacuum, and baryonic densities in units of the crit- 
ical density, (t| is the variance of linear fluctuations on 
the 8/i~^Mpc scale, and n is the "tilt" of the primordial 
power spectrum. The Eisenstein & Hu (1999) transfer 
function is used throughout. 

As in our earlier work (STDOl), simulations were con- 
ducted with a parallel OpenMP based implementation of 
the "HYDRA" code (Thacker & Couchman 2006) that 
uses the Adaptive Particle-Particle, Particle-Mesh algo- 
rithm (Couchman 1991) to calculate gravitational forces, 
and the SPH method (Lucy 1977; Gingold & Monaghan 
1977) to calculate gas forces. Gas densities and ener- 
gies are calculated using the standard SPH smoothing 
kernel method (for exact details see STDOl), with the 
kernel tuned to smooth over 52 particles; and radiative 
cooling is calculated using standard tables (Sutherland 
& Dopita 1993). We have kept the metallicity constant 
at Z = 0.05, to mimic a moderate level of enrichment in 
the galaxy formation process. However, this is an under- 
estimate of the intracluster metallicity, Z = 0.3, which 
seems to be an approximately universal value to inter- 
mediate redshifts (Tozzi et al. 2003). Lastly, because the 
epoch of reionization is poorly known, and because we are 
primarily focusing our attention on mass scales greater 
than 10^° Mq, we do not include a fiducial photoioniza- 
tion background in the simulation. 

We also do not include the so called "V/i" terms (Nel- 
son & Papaloizou 1994, Serna et al. 1996, Springel & 
Hernquist 2002) in our implementation of SPH. This is 
potentially a significant concern in this investigation as 
we will examine gas entropy. However, tests on an ex- 
panding spherical shell problem (see STDOl) using 1000 
particles and no artificial viscosity or cooling, show en- 
tropy conservation to be accurate at the 6% level, while 
the combined gravitational and hydrodynamic energy er- 
ror is around 1.5%. These findings are in broad agree- 
ment with the discussions presented in Hernquist (1993) 
and Springel & Hernquist (2002), where a small entropy 
conservation error was always accompanied by a larger 
energy error when integrating the evolution of the en- 
tropic function in the absence of V/i terms. While a 
lower error in the entropy is desirable, for the present 
phenomenological investigation and given the gains we 
get from using a code without Vh terms, we consider 
this error acceptable. 

We simulated a number of different box sizes and par- 
ticle numbers to quickly assess the accuracy and numer- 
ical resolution dependencies in our model. A single large 
simulation was then run for statistical purposes, allow- 
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ing us to probe the bright end of the luminosity function. 
The specifics of each simulation, including box size and 
resolution are given in Table 1. We note that attempt- 
ing to simulate the formation of the very brightest end of 
the luminosity function with sufficient resolution to track 
smaller mergers is a difficult task, due to the scarcity of 
these objects. This is the fundamental motivation behind 
our progressing to a simulation with 2 x 640'^ particles. 

2.2. Identification of Quasars 

A secondary motivation of this paper is to compare 
the simulation results directly with semi-analytic pre- 
dictions. We have therefore taken the outflow model 
of SO04 and adapted it to our simulation as closely as 
possible, although in some cases, which we highlight, it 
was either not possible for us to match this model ex- 
actly, or we have chosen to make well-motivated changes. 
For completeness, we reiterate the salient features of the 
SO04 model in our discussion below. 

While in our previous work it was sufficient to track 
group mass evolution to identify star forming regions, to 
evaluate the quasar luminosity function it is necessary to 
track mergers of groups. We have used the same method 
as STDOl for identifying groups, which relies upon the 
local baryonic density field to pinpoint centers of mass, 
and a spherical overdensity procedure applied to identify 
the baryon group. From the baryon group an estimation 
of the total halo mass is derived by multiplying by ^o/flb- 
The resulting mass distribution function is in close agree- 
ment with that derived from fricnds-of- friends, provided 
suitable limits for the baryon spherical overdensity are 
chosen. To track merger events we rely upon group la- 
beling, and we label a merger as having occurred when at 
least 30% of the accreted mass does not come from a sin- 
gle massive progenitor. This procedure means that the 
first groups to form are also treated as merger events. 
For each merger event that will be tagged as a quasar 
we store the details including position and redshift in an 
output file. 

Once a group has been identified as satisfying this cri- 
terion, the dynamical time associated with the cold gas 
disk which feeds the AGN and the mass associated with 
the black hole must be calculated. For a given redshift, 
z, and virial density Pv{z), the implied virial radius for 
a group of N gas particles with mass is 



4/3^p„(z) 
The circular velocity is then 



-T:Gp^{z)rl 



1/3 



1/2 



(1) 



(2) 



and the dynamical time, t^, associated with a cold disk 
of gas, is defined by 



td = 0.055 X ry/vc 



(3) 



Note that in order to maintain the same relationship be- 
tween outflow velocity and black hole mass as Wyithe 
& Loeb (2002) and SO04, we choose a slightly larger 
time than was used in these studies. This is because our 
numerical model also includes thermal energy input to 
establish the correct initial post-shock temperature ac- 
cording to eq. lO below. 



While the observed Mbh — relation (Merrit & Fer- 
rarese 2001; Tremaine et al. 2002) infers that the black 
hole mass scales as cr", where a ~ 4-4.5, the Mbh — Vc 
has a slightly steeper slope because the Vc — (Jc is shal- 
lower than linear (Ferrarese 2002). Thus we assume an 
Mbh — Vc relationship given by. 



bh = 2.8 X 10« (, 



, . (4) 

.300kms-iy ^ ' 

Lastly, the black hole is assumed to shine at its Edding- 
ton luminosity (1.2 x lO'^* ergs s~^ Mq~^) for the dynam- 
ical time, td- Note that these values are slightly different 
than in SO04, but the overall relationship between Vc 
and luminosity is the same. 

2.3. Outflow Implementation 

Each AGN in our simulation is assumed to channel a 
fixed fraction fraction, Ck, of its bolometric energy into 
a kinetic outflow. The amount of energy deposited into 
the outflow is then 



Ek ^ 1.2 X 10^^ efe 



M, 



hh 



Mr. 



td\ 

— ergs. 

s 



(5) 



As in SO04, we shall adopt — 0.05 throughout this 
investigation, which is consistent with other literature 
estimates [e.g. Furlanetto & Loeb 2001; Nath & Roy- 
chowdhury 2002). The majority of mass in the outflow 
at the resolution we can simulate will have come from 
material surrounding the cold gas group. Therefore, as 
in our previous work, we model the expanding outflow as 
a spherical shell outside of the virial radius of the system. 
While the assumption of a spherical shell is a signiflcant 
oversimplification, given the bipolar nature of outflows, 
it is worth recalling that within the intracluster medium 
in galaxy clusters a bipolar outflow will still launch an 
ellipsoidal cocoon of shocked intracluster gas (Begelman 
& Cioffi 1989). Hence, we place the expanding outflow 
at a radius 2rmr and re-arrange all gas below a density 
threshold of 2.5pvir within this radius, but outside ryir, 
into an expanding shell. The density threshold prevents 
us from redistributing cold gas, which is known to be 
very stable against incoming shocks in SPH simulations. 
Once we have established the amount of mass available 
to create the shell, Mg, the velocity of the shell, Vg, is 
calculated from, 



w, = 1.13 — Grrin———- > — 



1 



1/2 



M, 



'Ns 



(6) 



i=N+l 



where the second term on the right-hand side denotes the 
potential energy subtracted as particles are moved from 
their initial position to the shell. As in STDOl, we add 
an additional rotational velocity component to the shell 
so as to preserve the angular momentum. 

Note that the prefactor in eq. © is less than V2 as a 
fraction of is channeled into establishing the correct 
post-shock temperature in the outflowing gas. This heat- 
ing is particularly important since it will help determine 
the fraction of impacted gas that is able to cool within 
a Hubble time. Under the assumption that the shell be- 
haves like a strong shock, the postshock temperature, T^, 
is 



Ts 



mke (lkms-i)2 



(7) 
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TABLE 1 

Simulation parameters 



Run N 


Box Size 


Mass Resolution 


Softening Length 


Steps (z=l) 




(Mpc) 


(Mq) 


(kpc) 




0400 2 X 40^ 


35 h'^ 


1.2 X 10^0 


24 /i-i 


1340 


0800 2 X 80^ 


35 h'^ 


1.5 X 10^ 


12 h-^ 


3678 


1600 2 X 160^ 


35 h'^ 


1.9 X 10® 


6 


5420 


3200 2 X 320^ 


35 h-^ 


2.3 X lO'' 


3 /i-i 


7381 {z=2.5) 


1020 2 X 640^ 


146 fc-i 


2.2 X 10* 


9 h-i 


10420 {z=1.2) 



Whife the semi-anafytic modef in SO04 assumes that this 
heating appHes to the galaxy as well, here we only heat 
the material in the outflowing shell. Our motivation for 
this choice is the short cooling time of gas in galactic 
halos and, on a secondary level, the coUimated bipolar 
nature of the outflows. The radial expansion of the shell 
agrees with analytic predictions (STDOl). 

Since our resolution is insufficient to provide detailed 
knowledge of the inner structure of galaxies we imple- 
ment star formation on the basis of a merger model. 
Following a major merger we convert 10% of gas in the 
galaxy into stars particles. While this method is known 
to be a good model of high redshift star formation in low 
mass halos (STDOl) it does not track quiescent mode 
star formation which is the primary mode of star forma- 
tion in the higher mass galaxy population. We emphasize 
that the purpose of this study is to focus on the hydrody- 
namic evolution of the IGM and we are not attempting 
to calculate a luminosity function for galaxies, instead 
we use them largely as a tracer population. 

3. DISTRIBUTION AND EVOLUTION OF QUASARS 

3.1. Quasar Clustering 

In the following three subsections we study the opti- 
cal properties of quasars, as quantified in the rest-frame 
B-band. In keeping with our previous investigations, as 
well as the observations in Elvis et al. (1994), we relate 
the luminosity in this band to the overall bolometric lu- 
minosity by assuming a fixed ratio of LboI = 10.4Ls at 
all luminosities and redshifts. Fixing this value also al- 
lows for direct comparison with Wyithe & Loeb (2003). 
Finally, the LboI of the quasar associated with each out- 
flow is simply computed as LboI = Ek^ktd^ ■ 

We begin by addressing the spatial distribution of 
quasars, which serves as a check of our merger-based ap- 
proach. To quantify this distribution, we construct the 
spatial correlation function of quasars using the center- 
of-mass information from the outflow data produced in 
the simulation. In principle, this should be computed 
accounting for the finite lifetime of each quasar, ac- 
cording to eq. 0. In practice, these times are long 
enough that such effects can be ignored for distances 
< ctd ~ 20(1 + z)^^/^ comoving Mpc. Thus we cal- 
culate the 3-dimensional real space correlation function 
using the simplest estimator, 

g,,(z,m), + l= (8) 

where DD{z, m)k is the number of pairs with a magni- 
tude greater than some limit, separated by a comoving 



difference corresponding to a bin k, and RR{z, m)k is the 
average number of pairs that would be found at a given 
separation in a random distribution of points with an 
overall density equal to the mean density of observable 
quasars. 

We next adopt a fixed magnitude limit in the B-band 
of 20.84, to allow for comparisons with observations from 
the 2dF quasar redshift survey (Groom et al. 2001; 2002), 
which have an overall photometric b band limit of 20.9, 
where B f=i b + 0.06 (Goldschmidt & Miller 1998). This 
can be computed as 



B = 5.5- 2.51ogi 



is 



51ogi 



( d 



VlOpc 



+ 2.5(1 -a.)logio(l + z). 



(9) 



where d is the comoving distance to the quasar, and 
Ui, — —0.5 is the typical slope of the quasar power-law 
continuum (Wyithe & Loeb 2005). 

A detailed examination of short range quasar correla- 
tions is given in the Sloan binary quasar study of Hen- 
nawi et al. (2006). To compare to their sample we project 
the redshift space correlation function within a maximum 
velocity range of |Ai;| < 2000 km s^^, to give 



,(i?,z) 



^qq{R,S,z)ds. (10) 



Here the integration limits are set by Smax — —Smin — 
2000 /aH{z) corresponding to the width of the velocity 
interval. Note that in determining Wqq{R, z), rather than 
integrating over the redshift space correlation function 
we instead use the real space version. This is a good 
approximation since the defined velocity interval is suffi- 
ciently large to contain most of the velocity distribution 
in the redshift direction. 

In Figure ^ we plot the real space correlation func- 
tion in the z — 2.0 — 3.0 and z = 1.2 — 2.0 ranges. On 
scales larger than 1 Mpc the agreement with the Groom 
et al. (2001) results is strong, with our results showing 
a small excess over the mean observational result. The 
overlay of the higher redshift results in the middle panel 
shows, as expected, that the redshift evolution at these 
early times is very weak. On sub-Mpc scales there is 
a noticeable turn-up, which is best examined using the 
projected correlation function Wqq{R,z), shown in the 
bottom panel. To compare to the Groom et al. (2001) re- 
sults we have projected their real space correlation func- 
tion fit, and also extended this fit below their nominal 
cutoff. This provides a clear quantitative baseline for 
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Fig. 1. — Top: Real space correlation function of B < 20.9 
quasars from 2 = 2 to 3. The points are taken from our large 1020 
run, while the solid line is the (2.1 < 2 < 2.9) fit to the 2dF quasar 
redshift survey by Croom et al. (2001), bounded by the measure- 
ment errors. Note that this measurement does not extend to sep- 
arations below O.SMpc h~^. Center: Real space B < 20.9 quasar 
correlation function from the 1020 simulation from z = 1.2 to 2, 
as compared with the 1.35 < 2 < 1.7 Croom et al. (2001) results. 
Symbols are as in the upper panel, while the open circles are the 
simulated 2 < 2 < 3 correlation function, shown for comparison. 
Bottom: Angular correlation function of B < 20.9 quasars. The 
solid circles are taken from the simulation at 2: = 1.2 — 3.0, while 
the solid line is the Croom et al. (2001) observations, projected 
and averaged over the same redshift range. The cross-hatched re- 
gion is an extrapolation of the Croom et al. (2001) fit to smaller 
separations, which helps to highlight the excess in Wqq at small 
separations observed by Hennawi et al. (2005). This excess is well 
reproduced in our simulations, where it arises from gravitationally 
bound pairs of quasars (the so-called "one halo contribution"). 



examining the short-scale clustering excess observed by 
Hennawi et al. (2006). The precise position and magni- 
tude of the clustering excess are reproduced extremely 
well within our simulation, which we take as both sup- 
port for our approach and validation of the Hennawi et 
al. (2006) results. 

This turn up, which has also been observed in the 
high-redshift Lyman break galaxy population (Ouchi et 
al. 2005) as well as in a local sample of galaxies (Zehavi 
et al. 2004) , is most likely due to gravitationally bound 
pairs of quasars that are orbiting each other. This so- 
called "one halo" contribution {e.g. Bullock et aZ. 2002; 
van den Bosch et al. 2003; Magliocchetti & Porciani 2003) 
should become important at distances less than 



ilhalo 



(47r/3)r!oPc 



1/3 



(11) 



where mg is the mass of the halos associated with the 
quasars above our magnitude limit. From the large scale 
clustering this is w 2 x IO^^Mq, corresponding to the 
~ 1.0 Mpc position of the turn-up, lending further weight 
to this interpretation. 



3.2. Quasar- Galaxy Cross- Correlation Function 

As a complementary investigation, we also examine 
the cross-correlation function between quasars and galax- 
ies, S^gg. By cross-correlating these two populations we 
are directly able to evaluate whether galaxies containing 
quasars are clustered differently than similar-mass qui- 
escent galaxies. Early observational attempts to mea- 
sure ^qg were limited by sample size and thus a bias 
toward 2-dimensional angular measurements prevailed 
{e.g. Elhngson, Yee & Green 1991; Smith, Boyle & Mad- 
dox 1995; Croom & Shanks 1999). However, the SDSS 
and 2dF quasar surveys have enabled cross correlations 
in 3-dimensions below a redshift limit of z < 0.3 (Croom 
et al. 2003; Wake et al. 2004) and a study at intermedi- 
ate redshifts using the DEEP2 data has been undertaken 
(Coil et al. 2006). To date, these investigations have not 
uncovered any bias in £^gg on scales down to the minimum 
scale to which they are sensitive, which is around 1 Mpc. 

To evaluate ^gg we use the quasar catalog from the pre- 
vious section, combined with a galaxy catalog evaluated 
with a EOF group finder on the baryonic material in our 
simulation. We use a linking length b = 0.065 to find 
groups with an outer density limit oi S ~ 2000. With a 
baryonic mass cut of 10^°'^ Mq we find 37,995 groups, 
and for 1O^^M0 we find 11,857 groups. The estimator 
for the cross correlation function is 



iqg{z,m)k + 1 = 



DgDg{z,m)k 
RqRg{z,m)k ' 



(12) 



where DqDg{z, m)k is the number of quasar-galaxy pairs 
above the magnitude limit in bin k, and RqRg{z,m)k 
is the number of quasar-galaxy pairs that would be ex- 
pected if these objects were randomly distributed with 
the same densities as in our simulation. 

Our results are plotted in Figure 12 in which we now 
employ an absolute magnitude limit of Mb = —22, to 
better compare with observations. On scales larger than 
1 Mpc, S^gg is indistinguishable from ^^g, in agreement 
with observations. This is true regardless of whether we 
choose the M{, > 10^"-^ Mq or Mb > lO^M© galaxy 
populations. However, on scales below 600 h~^ kpc, £_qg 
exhibits a clustering enhancement. At first glance, this 
would appear to be consistent with our results from §3.1, 
which show that one-halo effects can create an excess of 
clustering at small scales. However, the explanation can- 
not be this straightforward. In the bottom panel of Fig- 
ure El we plot the ratio S,qg/£,gg, which shows explicitly 
the turn-over from large-scale agreement to short-scale 
excess. At large separations, active galaxies are clustered 
very similarly to the general population, consistent with 
the DEEP2 results. Yet at small separations, a dramatic 
change occurs. Despite the fact that the one-halo contri- 
bution is implicitly included in , the amplitude of the 
break in the cross correlation function exceeds that of the 
galaxy autocorrelation function by a factor of ~ 2.5. 

Thus it appears that our identification of quasars with 
mergers enhances their clustering on the smallest scales. 
At first it seems that this is strongly at odds with pre- 
vious theoretical studies of mergers (Percival et al. 2003; 
Scannapieco & Thacker 2003). However, these studies 
were targeted to separations larger than 1 Mpc, where 
the two-halo term is the dominant contribution to the 
correlation function. The excess we find here occurs 
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Fig. 2.— Top: Real space cross-correlation function, ^qg: 
of Mg < —22 quasars to galaxies with baryonic masses above 
IO^^'^Mq (circles) as compared to the autocorrelation function of 
these galaxies (squares). The selected mass and magnitude limits 
have been chosen to be in broad agreement with the quasars and 
galaxies used in measuring the cross-correlation function from the 
DEEP2 survey (Coil et a/. 2006). These observations are shown 
as the shaded region, which is an estimate that adopts the errors 
from the projected cross-correlation function, around the best fit 
to the spatial correlation function, ^qg(r) = (r/3.45)~^'®® . Finally 
the solid lines gives the linear correlation function of lO^'^ Mq ha- 
los. Center: Real space quasar-galaxy cross-correlation function 
and galaxy-galaxy correlation function, but now associating galax- 
ies with objects with baryonic masses above IO^^A/q. Symbols are 
as in the upper panel. Bottom: Ratio of quasar-galaxy cross- 
correlation function to galaxy-galaxy correlation function in the 
Mi, > 1O^°-^M0 case (squares) and Mf, > lO^^M© case (crosses). 
The cross correlation function not only exhibits a break below 1 
comoving Mpc , but this break is stronger than in the galaxy- 
galaxy autocorrelation function. The strengthening of this break 
is related to the merger nature of quasars in our simulation, and is 
consistent with local measurements (Serber et al. 2006), but unfor- 
tunately occurs at too small scales to be measured at 2 1 from 
current surveys. 



purefy in the one-halo regime, meaning that mergers have 
an excess of very close neighbors. 

In support of our results, local studies of quasar-galaxy 
clustering have found an excess of galaxies at radii < 0.5 
Mpc, completely analogous to the one in our simulation. 
Studying a, z < 0.3 sample of quasars drawn from the 
Sloan Digital Sky Survey (SDSS), Serber et al. (2006), 
have found that Mi < —23.3 quasars are more than 
three times more clustered than L* galaxies on <; 0.1 
Mpc scales, although they cluster similarly to L* 
galaxies on larger scales. While this study was carried 
out at much lower redshift, these quasars have roughly 
the same intrinsic magnitudes as those in Figure 2, as 
can be estimated assuming a typical redshift of 2; « 1.4, 
which for our simulated sample gives a distance modulus 



of w 45 or Mb < —24. Thus there seems to be mounting 
observational and theoretical evidence that the correla- 
tion function of the products of mergers, while only very 
weakly enhanced at large separations, may nevertheless 
be more strongly enhanced in the one-halo regime. 

In fact, an intriguing possibility is that this is caused 
by three-body interactions in which a third galaxy re- 
moves angular momentum from a nearby close pair. This 
suggests that dynamical friction may not always be the 
dominant process driving galaxy mergers. Rather, a sig- 
nificant number may be caused by a process more akin 
to the formation of tight binaries in dense star clusters 
{e.g. Rasio, Pfahl & Rappaport 2000). Clearly this issue 
merits future investigation. 

3.3. Optical Quasar Luminosity Function 

To construct the luminosity function for each luminos- 
ity and redshift bin, we calculate the number of quasars 
in this bin times the total time these objects are shining, 
and divide by the time interval, the width of the bin, and 
the volume of the simulation. That is for a given redshift 
bin i and a given luminosity bin j the luminosity function 
is simply 



'^'^^ ^ VAt.ALB. ^ 



td,k, 



(13) 



where the sum is over all quasars with redshifts and lu- 
minosities associated with the i, j bin, which spans a 
time interval Ati and a range of luminosities ALbj- In 
Figure|21we plot the resulting luminosity function for our 
fiducial 1020 simulation. 

The dotted line in this plot gives the Wyithe & Loeb 
(2003) estimate of the luminosity function, which is sim- 
ply based on a merger prescription and does not account 
for feedback. Comparing this estimate with our simu- 
lation results uncovers a clear turn-down in the number 
of Lb > IO^'^Lq quasars at z ^ 2. However, comparing 
the simulation results with the measured points make it 
immediately clear that this turn down is not as strong 
as seen in the observations. This means that the sup- 
pression is much weaker in the simulation than in the 
semi-analytic SO04 model, which was found to be a good 
fit to the observations. This is a surprising result which 
needs to be understood in more detail, and, to explore 
this further, we recalculate the luminosity function but 
imposing by hand the precise gas heating methodology 
used in the semi-analytic approach. 

As emphasized in Oh & Benson (2003) the ability of gas 
to cool is insightfully described from the perspective of 
entropy. Under the frequently used definition of entropy, 
S — T/n'^/^, the isobaric cooling time may be written, 



'^cool 



nlk{T) 



2/x2Ti/2A(r) 



(14) 

where /i — 0.62, /ig = 1.18, and F(T) serves as temper- 
ature dependent normalization. By equating the cooling 
time to the Hubble time, tjj, we can derive a critical en- 
tropy value: gas above this entropy limit is unable to cool 
within the Hubble time and is thus effectively removed 
from the galaxy formation process. The critical value is 
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Fig. 3. — Evolution of the B-band quasar luminosity function. Here the data points are taken from Pel (1995, open circles), which 
are derived from the compilation by Hartwick & Schade (1990), and Fan (2001, open triangles). Simulation results are given by the solid 
points, while the dotted line is the simple estimate from the analytic model by Wyithe & Loeb (2003), which does not include feedback. 
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raw simulation results, while the bottom row shows the luminosity function derived by applying the same physical model as SO04 to the 
simulation quasar catalog, along with the fiducial SO04 model with a 5% outflow efficiency (dashed line). The difference at low redshift 
between the two models arises because the semi-analytic model treats heating and cooling in a significantly different way to the simulation. 



(SO04), 

S^.it = (280 keVcni2)(l 4 

A(r„ 



E{z 



(1 + z)3/2 



2/3 







6.3 X 10-22 ergs s-1 cm3j ' ^^^^ 

where A(Tmin) corresponds to the minimum cooling time 
below 10^ K (Tmin ~ 2.3 x 10^ K), and E{z) = [f7„,(l + 

Having defined the concept of critical entropy, we first 
estimate the total amount of mass that can be heated 
above 5crit using eq. (17) of SO04, 

Afex(<5, Z, Mbh) = 4.6 X 10^2 Mq 5i00,crit(z)"' 

xEeoS-^^Ul + zr\ (16) 

We then search to a radius R — min{Rs, Rheat), where 
Rs is the shell radius versus time given by, 

Rs = 1.7MpcEli' d-'/' (1 + z)-3/5 t'J^^, (17) 

and Rheat is the radius of the region heated above 5'crit 
as calculated in the SO04 model 

Rheat = 5.6Mpc5ioo,crit(^)~'/'-B6o'C'^' (1 + ^)"'^' ■ 

(18) 



Our assumed ^dit is 60 keV cm^, corresponding to the 
metallicity value of Z = 0.05 used in our simulation, and 
we taken an average post-shock overdensity, Ss, of 20. If 
a quasar is found within R, we subtract Mex from the gas 
mass associated with it, and if Mex exceeds this mass, we 
remove it from our catalog altogether. 

The resulting luminosity function is plotted in the 
second row of Figure |21 In this case, there is sig- 
nificantly better agreement between the data (and the 
semi-analytic model) for this revised luminosity function. 
While this model does not precisely reproduce the strong 
knee observed in the SO04 model and the magnitude 
of the turn-down at higher luminosities, the improved 
agreement is compelling. It is thus clear that the pri- 
mary cause of the difference between these models, is 
not their varying approaches to modeling quasars them- 
selves, but rather in the simplified exclusion conditions 
as described by eqs. ((TCll - ((T^ . 

While these equations capture most of the salient fea- 
tures of shock heating, there are two major effects that 
they fail to address. Firstly, they do not differentiate 
between material within galaxies and material in the in- 
tergalactic medium, while in our simulation we do not 
apply the outflow heating to the host galaxy itself. Sec- 
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ondly, eqs. H16|) - (|18|l assign a single density to aU the 
gas associated with an object that is overtaken by an 
outflow, whereas heating processes in the simulation are 
directly affected by the detailed substructure in this gas. 

To explore the relative impact of these two effects, we 
examine the distribution of progenitor halo labels in the 
three largest outflows in our simulation at an epoch of 
z ~ 1.3. These systems have baryonic masses ranging 
from 4.7 x lO^^ Mq to 5.1 x lO^^ Mq, and the majority 
of particles within them have a single label associated 
with the previous outflow event. To determine the num- 
ber of particles with dissimilar labels we place a sphere 
at the center of mass of the group and then fit a radius 
(by hand) to enclose the outer boundary of the particles 
with the main group label. The typical radius of this 
outer boundary was 250 kpc. For all systems we found 
a significant amount of substructure as evidenced by all 
groups having at least 18, and typically more than 40, 
distinct labels each with at least 20 particles. The total 
number of particles with labels distinct from the main 
group was always close to the merger mass limit, indica- 
tive of a system just about to undergo the outflow event 
associated with the merger. 

This results seems to suggest that both effects may 
be contributing to the reduced suppression, as each out- 
flow is associated with a large group with a single index 
(that might have been disrupted if in-galaxy heating was 
included), which merges with a collection of remnants 
with many indices (that might not have been accreted 
if in-shock heating were more efficient). Ultimately, the 
stronger agreement for the revised catalog versus that of 
the original simulation indicates that both of these is- 
sues need to be explored further. On the simulation side 
this would involve contrasting the present results with a 
model that includes more efficient ejection of gas from 
galaxies themselves. Likewise, in the semi-analytics the 
efficiency of heating could be parameterized by adding 
an additional parameter to account for in-shock cooling, 
although it would be necessary to conduct detailed simu- 
lations of this process to precisely calibrate this number. 
These are significant issues which we will return to in the 
discussion. 

Finally, as a test of convergence, in Figure 0] we com- 
pare the luminosity function in our fiducial simulation 
with that derived from our 35 Mpc^ simulations. In 
this case we do not impose eqs. (|16|l - (|18ll . At all red- 
shifts we obtain good agreement between runs over the 
luminosity range spanned by these smaller simulations 
volumes, although the 400 results are very noisy. Note 
that while in the z k, 2.25 column, the luminosity func- 
tion in the 3200 run has fewer large quasars than the 
other runs, this is due to small number statistics and the 
fact that it was stopped earlier than the other simula- 
tions due to the excessively large number of time-steps 
required. Similarly, Lb > IO^'^'^Lq quasars are so rare 
that they can not be compared between the 1020 simu- 
lation and the test simulations, motivating our choice of 
an extremely large volume for this run. 

3.4. Hard X-ray AGN Luminosity Function 

Our simulation can also be directly compared with ob- 
servations of the hard X-ray luminosity function, which 
we construct from our bolometric luminosities using a 
similar approach to that for the optical quasar luminosity 



function. To convert from the bolometric luminosity to 
Lx we use results from Marconi et al. (2004) who calcu- 
lated the expected correction for k, 10^^ Lq AGN with a 
spectral template motivated by recent observations. The 
ratio of the bolometric luminosity to the hard X-ray band 
(2-10 keV) is given by a third-degree polynomial, 

log[L/L(2-10keV)] = 1.54+0.24/:+0.012/:2+0.0015/:^ 

(19) 

where L = log(iBoi) — 12, and XboI is given in Lq. 

Early observational work used ASCA data (Boyle et 
al. 1998) and BeppoSax data (La Franca et al. 2002) to 
show that the hard X-ray luminosity function (HXLF) 
was evolving strongly between z = Q and 1.5, consistent 
with pure luminosity evolution. More recently, Cowie et 
al. (2003) used Chandra data in two redshift bins (2=0.1- 
1 and z=2-4) to argue that the AGN number density for 
luminosities lower than lO*"' erg s~^ seems to peak at a 
lower redshift than those of higher luminosity. This an- 
tiheirarchical evolution was demonstrated definitively by 
Ueda et al. (2003), who carried out a comprehensive com- 
pilation of HEAOl (Piccinotti et al. 1982; Grossan 1992), 
ASCA (Ueda 2001, Akiyama et aZ. 2003) and Chandra 
(Brandt et al. 2001) data to derive the comoving spa- 
tial density of AGNs in three luminosity ranges between 
log(Lx) = 41.5 and log(Lx) = 48. 

In Figure [S] we compare our derived spatial densities 
to the observational data in the \og{Lx) = 43 — 44.5 and 
\og{Lx) = 44.5 — 48 luminosity bands in this sample. 
Both the qualitative and quantitative predictions of the 
simulation agree with the measurements: the downsiz- 
ing trend is apparent in both luminosity bands and the 
overall normalizations agree well. However, below z ~ 2, 
we are faced with two minor issues. Firstly, in the lower 
luminosity bin it appears that the luminosity function 
is turning down slightly too quickly, as the observations 
suggest that the turn down in this luminosity range oc- 
curs after z — 1. Secondly, the brightest band, while 
turning down at the observed epoch of z = 2, does not 
perfectly follow the observational trend. Imposing the 
exclusion conditions represented by eqs. ((TE|l - lfTF|l im- 
proves this fit somewhat, although these differences are 
small compared to the measurement errors from the ob- 
servations. Imposing these conditions has no impact on 
the lower luminosity bin. In general, these results are 
consistent with those in H3.3I 

4. OTHER IMPLICATIONS 
4.1. Impact on the Intergalactic Medium 

While our study has been focused on the properties 
of the quasar population itself, our simulations naturally 
have predictions for the more tenuous gas surrounding 
large galaxies. In fact, as discussed above, the most clear 
observational evidence for widespread nongravitational 
heating lies not in the galaxy population, but rather in 
the properties of the diffuse gas in galaxy clusters. 

To examine the impact of outflows on this material 
we first focus on the total amount of gas that has been 
shocked to 5 > 5crit, such that it no longer participates 
in fueling further generations of quasars. The redshift 
evolution of the mass fraction of this gas is plotted in 
the upper panel of Figure El which shows that quasar 
feedback is primarily a low- redshift phenomenon. Thus, 
above 2^3 less than 3% of the gas in the simulation has 
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Fig. 4. — Evolution of the B-band quasar luminosity function as a function of resolution. Again the data points are taken from Pei 
(1995, open circles) and Fan (2001, open triangles), as compiled by Hartwick & Schade (1990). Simulation results are given by the solid 
points, while the dotted line is the simple estimate from the analytic model by Wyithe & Loeb (2003), which does not include feedback. 
From left to right the columns give results at redshifts of 1.0-1.75 (1.2 - 1.75 in run 1020), 1.75-2.75 (2.5-2.75 in run 3200), 2.75-3.75, and 
3.75-4.75 respectively. From top to bottom, the rows show results from our large (146/i~^ Mpc^) run 1020, and the smaller (35?t~^ Mpc^) 
runs 0400, 0800, 1600, and 3200. Symbols are as in Figure ISl and, for comparison, the crosses in 2 = 2.5 — 2.75 panel of the 3200 run show 
the results of the 1600 run limited to this same redshift range. 



been affected, consistent with the lack of suppression of 
the luminosity function at these redshifts. At lower red- 
shifts, however, the S > ^crit rnass fraction grows prodi- 
giously, preventing roughly 20% of the gas in the simula- 
tion from cooling. This is consistent with the turn-down 



in the luminosity function seen in Figure |31 which while 
not as efficiently quenched as the semi-analytic results, 
nevertheless differs substantially from the pure-merger 
predictions. 

As a test of convergence, we also plot in this panel the 
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redshift redshift 



Fig. 5. — Hard X-ray luminosity function. The solid and 
dashed lines give the number density of X-ray luminous (44.5 < 
logio(Lx) < 48) and X-ray faint (43 < logio{Lx) < 44.5) AGN in 
our 1020 simulation, respectively. In the luminous case, the lower 
solid curve shows the results in which the exclusion conditions, eqs. 
I16i - 1181 are imposed. These models are compared with obser- 
vations of luminous (circles) and faint (squares) AGN as compiled 
by Ueda et al. (2003). 



S > S'crit mass fraction from each of our smaller simula- 
tions. These range from the 0400 run, in which particles 
are 64 times more massive than in the 1020 simulation, 
to the 1600 run, in which particles are 0.125 times the 
mass of those in the 1020 run. As increasing resolution 
adds a large number of low-mass, high-redshift outflows, 
the mass fractions at high redshift increase monotoni- 
cally with resolution. At lower redshift, however, the 
mass fractions approach each other asymptotically, and 
in the important z <^ 3 range, this quantity is largely con- 
sistent across runs. However, this mass convergence does 
not give a complete picture of the effect of resolution. 

In the lower panel of Figure 1^1 we plot the evolution of 
the volume filling factor of S > ^crit gas in each of these 
runs. To calculate these quantities in the SPH method 
it is necessary to first smooth the particle data on to a 
grid. In the 1020 case we do so on a 1340"^ mesh, so that 
the smoothing scale for the filling factor at expansion 
factor a is 0.155a Mpc, which is considerably above our 
minimum smoothing length. In the other runs we use a 
mesh of size twice the particle resolution, for example, 
the 2 X 160^ run was smoothed on to a 320'^ mesh. 

The data point in this panel gives the upper limit of on 
the volume filling factor provided by the Lyman-a forest 
(Levine & Gnedin 2005), which is well above our results, 
as expected from the semi-analytic estimates in SO04. 
For comparison, we also plot the results from a range 
of models taken from the N-body + semianalytic study 
of Levine & Gnedin (2005). While this is substantially 
higher than our results, this is to some degree due to 
the fact that they computed the full volume impacted 



Fig. 6. — Top: Mass fraction of gas heated above Scrit as a 
function of redshift. The solid line shows the results from our large 
1020 run, while the short dashed, long dashed, dot short-dashed, 
and dot long-dashed, lines (which move from higher to lower filling 
factors) show the results of our 3200, 1600, 0800, and 0400 runs, 
respectively. Bottom: Volume filling factor of gas heated above 
Scrit as a function of redshift. Lines are as in the upper panel, 
while the point gives the hya forest constraint discussed in Levine 
& Gnedin (2005, see also Dave et al. 1999). Finally, the shaded 
region gives the estimate of the volume filling factor impacted by 
AGN outfiows in Levine and Gnedin (2005), and is bounded from 
below by their tk = 0.05, tagn = 10^ yr model and bounded from 
above by their ej, = 0.05, tagn = 10* yr model. 



by quasars outflows, rather than only the volume heated 

above S'crit- 

Furthermore, it is clear from this figure that the vol- 
ume filling factors are significantly different across simu- 
lations, even at the lowest redshifts. This result seems at 
odds with our mass-fraction measurements. To explore 
this issue further, in Figure[3we choose a fixed redshift of 
z = 3 and plot the mass fraction and volume filling frac- 
tion above a threshold entropy Sthrcsh, which we allow 
to vary. For all three models, the mass fraction is only a 
weak function of Sthrosh for all entropy values near 5crit- 
This means that small differences in entropy have only a 
small impact on the number of particles prevented from 
cooling, and therefore both the S > S'crit mass fraction 
(shown in Figure and the luminosity function (shown 
in Figure^ are similar across runs. Essentially, at z = 3, 
the particles are divided into two types, those whose en- 
tropies are well above Scrit , and those that are far below 
this critical value. 

In the lower panel of Figure 0] we plot the volume fill- 
ing factor as a function of Sthresh, again at z = 3. In 
this case, near Scrit the volume filling factor is a strong 
function of our threshold entropy. This suggests that the 
high entropy gas is largely found in low density environ- 
ments, so that changes in S around Scrit pass through a 
region where the volume can change rapidly, but there is 
actually little mass to modify the overall mass fraction. 
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Fig. 7. — Top: Mass fraction of gas heated above a given 
entropy threshold, S'thrcsh in or 3200 (short-dashed), 1600 (long- 
dashed), and 0800 (dot-dashed) simulations, at a fixed redshift of 
2 = 3. Center: Volume filling factor above a given threshold at 
2 = 3. Lines as in the upper panel. For our choice of metallieity, 
'S'crit(^ = 3) = 19 keV cm'^, and near this critical value, the vol- 
ume filling factor is a strong function of the threshold entropy, but 
the mass fraction is almost constant. Bottom: Average differential 
overdensity (change in mass fraction over change in volume filling 
factor) as a function of 5thrcsh- 



In the lower panel of Figure 01 we plot the the differ- 
ential overdensity, AM/AV, that is the change in the 
mass fraction over the change in the volume filling fac- 
tor, as a function of S. This confirms that the majority of 
S « ^crit gas is in environments only a few times denser 
than the mean, and that the density of this material is 
increasing strongly as a function of entropy. 

FigurelHlillustrates why this is this case. Here we show 
the entropy distribution in slices taken from our 1020 
simulation at the final output redshift. It is clear from 
this plot that the boundary between S > ScHt and sub- 
critical gas (indicated by the white lines) lies at the very 
edges of cosmological halos, where the density is drop- 
ping off strongly. These boundaries are defined by com- 
paratively few SPH particles, and thus their positions 
can change rapidly following small changes in the par- 
ticle distribution. Finally, the thin slice shown in the 
smallest-scale panel in Figure |S1 uncovers the presence 
of 5 < Sent subclumps within a larger (high entropy) 
region. Again, the presence of this cooling substruc- 
ture was not included in the SO04 models, and is one 
of the key causes of the differences between these semi- 
analytical results and those presented here. 

4.2. Impact on Clusters and Groups 

Recent Chandra observations of Hydra A (Nulsen et al 
2005) have uncovered evidence for AGN activity within 
clusters at intermediate redshifts. While cooling flow 
clusters at lower redshift seem inconsistent with the idea 



of powerful shocks {e.g. Voit & Donahue 2005, Croton et 
a/. 2006), it remains an intriguing possibility that these 
systems underwent an earlier period of strong feedback 
and have now settled into a quiescent state. The observa- 
tions of radio quiet clusters by Donahue et al. (2005) that 
show extremely long cooling times despite an absence of 
inferred black hole activity, are broadly consistent with 
this hypothesis. These results prompted the numerical 
investigations of Sijacki & Springel (2006) who showed 
that their model of AGN activity has comparatively little 
effect on the cluster Lx — T relationship when all mate- 
rial within the virial radius is included. Therefore in this 
section we examine the effect of our outflow model on 
the Lx — T relationship and the cluster entropy profile. 
It is worth noting that as our raw luminosity function 
shows we have an excess of bright quasars at z=1.25, the 
results we derive in this section are an upper bound on 
the effect of outflows on clusters. 

To find cluster groups we first begin from a friends- 
of-friends b=0.2 catalogue of the 1020 simulation. The 
centers of mass are evaluated for these groups, and then 
used as the beginning stage of an iterative spherical- 
overdensity group finder that searches radially outward 
until the group is below the density threshold. The center 
of mass is then evaluated and used as the beginning point 
for the radial search, with the process being repeated 
five times. This technique has the advantage of bias- 
ing against mergers since the center of mass of mergers 
is usually sufficiently offset from the merging groups to 
stop the spherical-overdensity convergence process early, 
and the group is then discarded due to the low amount 
of mass found. We find 1272 groups by this process. 

To estimate the bolomctric luminosity of the simulated 
clusters we use 



No 



(piTLp) 



(20) 



and the emission weighted temperature is given by, 
_J2,m,p,A{T,)Ti 



(21) 



In these formulae, A(Ti) is the pure bremsstralung esti- 
mator of Pearce et al. (2000) (see also Navarro et al. 1995 
and Muanwong et aL2001, for similar approaches), and 
mi,pi,Ti are the mass, density and temperature of par- 
ticle i. While we could use the emission curve associated 
with the Z — 0.05 metallieity gas used in the simulation, 
the strong peak caused by coUisional excitation of He"*" 
at 10^ K will give a very high weighting to the cores of 
clusters, which in fact would probably have cooled sig- 
nificantly if we were using a, Z = 0.3 metallieity gas. 
Overall, as emphasized in Muanwong et al. (2001), our 
choice of a pure bremsstralung estimator will bias our 
luminosities low. 

Under the assumption of self-similar evolution 
for spherically symmetric clusters, and also that 
bremsstralung emission dominates in the X-ray band, the 
luminosity-temperature relationship scales with redshift 
according to (Kaiser 1986, Maughan et aZ. 2006), 

E{zy^A{z)-^^^L cx kT^fg, (22) 

where E{z) is the cosmological evolution factor, A corre- 
sponds to the virial overdensity at a given redshift {e.g. 
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Fig. 8. — Top: Contours of entropy from z = 1.2 slices through in our 1020 simulation. In the upper two panels the thickness of the 
slice is 2.4 Mpc h~^, while, to emphasize substructure, the thickness in the 36 Mpc and 18 Mpc panels is taken to be 1.2 Mpc 
and 0.6 Mpc h~^, respectively. All units are comoving. In each panel the white lines demarcate the boundary at which S = S^rit- In 
general, this occurs at the very edge of halos, where the density is dropping rapidly, but it also occurs in isolated subclumps in otherwise 
S > 5crit regions. 



see Bryan & Norman 1998) and fg is the chister gas frac- 
tion (which is assumed to be independent of z and kT 
following observations of high redshift clusters, e.g. Allen 
et al. 2002). Thus the redshift evolution in the Lx - T 
relationship can be scaled by dividing the luminosity by 
£;(z)(A(z)/A(0))-i/2. While we could thus scale the 
Lx — T relationship for z — clusters back to our final 
redshift, a recent analysis by Maughan et al. (2006) of 
the WARPS clusters in the region 0.6 < z < 1.0, pro- 
vides an unsealed Lx — T relationship at z « 1. Thus 
in what follows we use their results as a comparison. We 
also note that there are no observations of galaxy groups 
at these epochs. 
We plot our results for the Lx — T relationship in Fig- 



ure 13 and also show the fit of Maughan et al. (2006). It 
is immediately noticeable that our raw data show a very 
large scatter in luminosities. Since the luminosity of clus- 
ters is linearly weighted by the density, we decided to plot 
radial profiles of our clusters to determine whether any 
"overcooling" effects {e.g. Thacker et aL 2000) might be 
present (see figure [TT|l . The radial profile shows that the 
most massive clusters do indeed show a sudden upturn 
in density in their cores along the lines observed in test 
problems. We therefore applied a density-cut (filter) at 
^ < 5 X 10'' to the gas in our clusters, which cuts out 
this problem region. The resulting data is plotted in the 
right hand panel of figure El and shows a much smaller 
scatter. A least squares fit for our 1 keV and brighter 
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kT / keV 

ew / 

Fig. 9. — Lx — T relationship for the simulated clusters. On the left is the raw unfiltered data including the overly bright core regions 
(with overdensities exceeding 5 X 10^), while the right panel is the filtered data with these regions removed. Dotted lines correspond to 
least squares fits for clusters above 1 keV. The shaded area corresponds to Maughan et al. (2006) fit to the WARPS sample, with the upper 
limit set by taking the maximum normalization, and the lower limit by matching the exponent of our best-fit to the filtered data. 



unfiltered cluster catalogue gives 

/ rp \ 3.58 

Lboi = 2.82 X 10^2 (■ ^ j ^ (23) 
while the filtered catalogue gives, 

/ rp X3.18 

LboI = 2.09 X 10^2 (■ ^ j ^ (24) 

which is closer to, albeit slightly less luminous than, the 
Maughan et al. (2006) best fit of 5.4 x lO^'^iTen, /^e^f-^'^ ■ 
Overall, these results are strongly supportive of the con- 
clusions of Sijacki & Springel (2006). 

Lower redshift observations of groups [e.g. Xue & Wu 
2000), show the power law exponent for groups is close 
to 5. For systems below 1 keV (regardless of whether or 
not they are core filtered) we do not observe any steep- 
ening of the Lx — T relationship, and there is evidence 
of the relationship becoming shallower, indicative of the 
gas in this systems being quite strongly perturbed. Thus, 
rather than analyzing relaxed systems, we an are in fact 
analyzing groups that are radiating significantly due to 
the presence of an outgoing shock. In the event that 
this shock is sufficient to heat gas above ^crit, we might 
well expect, in the absence of significant further accre- 
tion and AGN activity, that these groups expand and 
cool before z — Q. To quantify this hypothesis, we plot 
the cluster and core entropy versus temperature in Figure 
1101 Our results are in broad agreement with the analy- 
sis of nearby clusters (Ponman et al. 2003, Finoguenov et 
al. 2005), and do indeed show that groups tend to show 
an excess of entropy in their cores. Therefore, we tenta- 
tively suggest that the X-ray emission of galaxy groups 
at 2; > 1 may well be a "smoking gun" for the AGN 
heating hypothesis. 




0.4 0.6 0.8 1 2 4 6 



kT,. / keV 

Fig. 10. — Entropy versus emission weighted temperature for 
different radial cuts. The 4-pointed crosses correspond to the en- 
tropy within the virial radius and are compared to the {z 0) 
self-similar (5 oc T) fit of Ponman et al. (2003) for the entropy 
within the isoo region, with the normalization taken to match the 
eight hottest clusters in their sample. The solid squares denote the 
entropy within 0.1r200i a measure of the cluster core entropy. The 
dashed line is the Finoguenov et al. (2005) weighted orthogonal re- 
gression fit to the Ponman et al. (2003) results. Our results show a 
significant turn-up in the core entropy at groups scales (< 1 keV). 



As well as the radial density profiles, we also show the 
temperature and entropy profiles in Figure [TTl The only 
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Fig. 11. — Density, temperature and entropy profiles for the four largest relaxed clusters at z = 1.2 in physical Mpc. A 100 particle 
moving average has been used to smooth the data. 



unusual feature in the temperature profiles is that the 

second ckistcr from the left has a very slightly inverted 
temperature profile following an extremely strong out- 
flow event. The resulting entropy profile is S" cx r^'"^, 
while the remaining profiles all match the 5* cx r^-^ pro- 
file reported elsewhere {e.g. Voit, Bryan & Kay 2005). 

5. DISCUSSION & CONCLUSIONS 

We have presented results from a suite of cosmologi- 
cal simulations that self-consistently follow the evolution 
of quasars and the outfiows associated with them. By 
tracking the merger history of halos and applying the 
quasar model of Wyithc & Locb (2003), we have been 
able to make direct predictions for the spatial distribu- 
tion and luminosity function of these objects. Our results 
are in excellent agreement with the observed correlation 
function of quasars on both small and large scales, re- 
producing both the power-law behavior measured by the 
2DF redshift survey on > 1 Mpc scales, and the strong 
break measured from the SDSS on < 1 Mpc scales. 

Furthermore, we predict that the quasar-galaxy cross- 
correlation function should show a small scale up-turn 
relative to the galaxy-galaxy correlation function. This 
occurs on sub Mpc scales, and therefore can only be 
measured from a very large sample of quasars, making 
its detection difficult even in the large DEEP2 survey. 
However, it is worth noting that the Dark Energy Sur- 
vey will produce a quasar catalog spanning 5,000 square 
degrees which combined with photometric redshifts for 
300 million galaxies should be able to suppress statis- 



tical uncertainties to a sufficiently low level to measure 
this upturn. We also note that this clustering excess is 
in qualitative agreement with the increase in the inte- 
grated galaxy overdensity at small scales for the SDSS 
quasar sample (Serber et aZ. 2006). Although the origin 
of the excess is uncertain, there is an interesting possibil- 
ity that the presence of an ancillary galaxy can accelerate 
the merger process via a three-body interaction. Given 
the extended nature of the mass distribution associated 
with galaxies it is not clear whether this mechanism will 
work in a similar way to stellar interactions. We plan to 
investigate this intriguing idea in the near future. 

As the correlation function is dominated by more pop- 
ulous low-luminosity quasars, our spatial results should 
be largely interpreted as lending support to our dark- 
matter modeling choice of merger model. Similarly, as it 
tracks the number of quasars formed as a function of the 
gas mass accreted onto galaxies, the luminosity function 
is more sensitive to the details of our feedback modeling. 
In this case, our simulation qualitatively reproduces the 
observed anti-hierarchical behavior, but the turn-down is 
much weaker than observed. Matching the suppression at 
the brighter end requires that we increase the efficiency 
of heating from the AGN outflows to mimic that assumed 
in semi-analytic models. This is equivalent to suppress- 
ing one key physical process, namely in-shock cooling in 
the presence of substructure, and including the ejection 
of gas from quasar host galaxies. These results empha- 
size how sensitive the luminosity function is to issues in 
the baryon physics and how the treatment of these issues 
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in semi-analytic models is still quite approximate. 

Investigation of the mass fractions and filling factors 
of gas above S'crit showed good convergence in the mass 
fractions at low redshift, but less so in the filling fac- 
tors. The differences between runs are most noticeable at 
high redshift, where successively higher resolution leads 
to the first generation of AGN outflows occurring at ear- 
lier epochs. While mass fractions show fairly strong con- 
vergence below 2 ~ 2, the tendency of gas above iScrit 

to 

occupy low overdensity regions makes an accurate calcu- 
lation of the volume filling factor difficult due to sam- 
pling issues. Those results are also further complicated 
by the known dependence of shock resolution on particle 
number {e.g. Thacker et al. 2000) where accurate model- 
ing of shock jumps in spherical collapse is reached once 
Ncoiiapse > 30,000. Nonetheless, the results do elucidate 
that the impact of quasar outflows is largely felt at low 
redshifts, in agreement with the downsizing trend. Per- 
haps the most interesting issue we have not yet explored 
with regards to filling factors is the relative impact of 
including a more bi-polar outfiow. However, there is no 
reason to expect a difference beyond a factor of two as the 
gas will flow toward low overdensity regions, as observed 
in supernova outflow calculations (STDOl). 

In somewhat denser environments, our results for the 
cluster Lx — T relation are in broad agreement with ob- 
servations. This is especially encouraging as our cluster 
modeling is not as sophisticated as other more targeted 
studies. In this case, the predominant view is that AGN 
heating is best modeled in terms of "hot bubble" ejec- 
tion from the brightest cluster galaxy (BCG) which then 
mixes with the surrounding material. Interestingly, the 
exact nature of this mixing is not well understood, and 
substantial differences are fomid between Eulerian and 
Lagrangian simulations, which by definition differ sig- 
nificantly in their treatment of advection. While SPH 
simulations inhibit the development of instabilities that 
promote mixing due to the necessary use of flow interpen- 
etration suppression, it has the advantage of exhibiting 
significantly less numerical diffusion than many Eulerian 
methods. Ultimately, input from laboratory experiments 
may well be necessary to help determine the correct mix- 
ing behavior. We also note that due to our inefficient star 
formation model, our simulations include a significant 
amount of cold gas in BCGs that will tend to promote 
an "overcooling" instability, which others have avoided 
by applying phase decoupling (Pearce et al. 2000) . We 
believe that this is a significant contributor to the lack of 
an entropy core in our radial profiles, and removing this 
central core produces an Lx — T relationship that is sig- 
nificantly less noisy while not exhibiting a large change 
in normalization at intermediate cluster mass scales. Ad- 
ditionally, the lack of a strong turn-down in the quasar 
luminosity function at z = 1.2 also promotes extremely 
strong outfiows in these clusters that may well be trans- 
porting high entropy gas from the core out to the edges 
of the cluster without significantly raising the entropy of 
the gas immediately surrounding the BCG. 
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Perhaps the most intriguing result to come out of our 
Lx — T study is the prediction that the gas in z w 1 
galaxy groups should be strongly perturbed by AGN ac- 
tivity, which is in the process of turning-off at that epoch. 
This effect is confined to small and high-redshift clusters 
for two main reasons. At lower redshifts, the comparative 
paucity of AGN activity will allow group gas to evolve 
largely adiabatically, decreasing Lx dramatically to es- 
tablish the steep Lx — T relation {Lx oc T^) observed 
locally. In more massive 2: « 1 clusters, on the other 
hand, the effect of outflows is largely masked when aver- 
aging over the material within r2oo- Thus galaxy groups 
at z « 1 represent the key mass and redshift range at 
which the AGN heating hypothesis is most likely to be 
testable though observations. 

Lastly, our simulation results indicate that simplifica- 
tions in current semi-analytic models may well be down- 
playing critical physics. While much attention has been 
paid to drawing broad conclusions from comparisons 
between observations and these models, the differences 
we have uncovered in this investigation are troubling. 
While, as expected, post processing of our simulation re- 
sults was able to reproduce the semi-analytic behavior, 
it is clear that in-shock cooling due to substructure in 
low-overdensity halos is an issue that must be consid- 
ered carefully in semi-analytic calculations. Fortunately, 
constraining the effect by simulation would not be diffi- 
cult, and ultimately a single parameter could be used to 
quantify this behavior. On the simulation side, it also is 
clear that a more detailed understanding of gas ejection 
from quasar hosts will be necessary for definitive conclu- 
sions to be reached. Nonetheless, the insight gained from 
this initial simulation study clearly serves to highlight 
the promise of a self-regulating picture of quasar forma- 
tion. The simple merger model of Wyithe & Loeb (2003), 
when supplemented with an outflow model, appears to 
represent a significant step forward in understanding the 
evolution of quasar clustering and the cause of their an- 
tiheirarchical low-redshift turn-off. 
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